Nonlinear elastic stress response in granular packings 
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We study the nonlinear elastic response of a two-dimensional material to a localized boundary 
force, with the particular goal of understanding the differences observed between isotropic granular 
materials and those with hexagonal anisotropy. Corrections to the classical Boussinesq result for 
the stresses in an infinite half-space of a linear, isotropic material are developed in a power series 
in inverse distance from the point of application of the force. The breakdown of continuum theory 
on scales of order of the grain size is modeled with phenomenological parameters characterizing the 
strengths of induced multipoles near the point of application of the external force. We find that 
the data of Geng et al. [l[ on isotropic and hexagonal packings of photoelastic grains can be fit 
within this framework. Fitting the hexagonal packings requires a choice of elastic coefficients with 
hexagonal anisotropy stronger than that of a simple ball and spring model. For both the isotropic 
and hexagonal cases, induced dipole and quadrupole terms produce propagation of stresses away 
from the vertical direction over short distances. The scale over which such propagation occurs is 
significantly enhanced by the nonlinearities that generate hexagonal anisotropy. 

PACS numbers: 45.70.Cc, 62.20.Dc, 83.80.Fg 
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I. INTRODUCTION 



The response of a granular medium to a localized 
boundary force has been investigated both experimen- 
tallyand numerically i, i, [li, i 0, i, U M M, 
E M M M [iGj . Experiments have shown that in 
disordered packings stress response profiles consist of a 
single peak that broadens linearly with depth 0, Q . For 
hexagonal packings of disks [l|, [2| or face-centered cubic 
packings of spheres 0, @1 j on the other hand, the stress 
response develops multiple peaks that seem to coincide 
with propagation along lattice directions. In two dimen- 
sions, a hexagonal packing is indistinguihable from an 
isotropic one in the context of classical (linear) elasticity 
theory [13, [3- Thus the observation of response pro- 
files in two-dimensional disordered and hexagonal pack- 
ings that differ significantly on scales up to 30 grain 
diameters (T], [2| requires consideration of nonlinear ef- 
fects. More generally, the applicability of classical elas- 
ticity togranular media is a question of ongoing research 

Classical elasticity for an isotropic medium predicts a 
single-peaked pressure profile that broadens linearly with 
depth Numerical results (see Ref. [13, for exam- 

ple) demonstrate responses well described by this solu- 
tion in regions far from a localized force in the bulk of 
a disordered frictional packing with more than the criti- 
cal number of contacts required for rigidity (the isostatic 
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point). Recent work by Wyart jl9[ and EUenbroek 
clarifies the onset of elastic behavior as average coor- 
dination number is increased above the isostatic limit. 



* Current address: Instituut- Lorentz, Universiteit Leiden, Postbus 
9506, 2300 RA Leiden, The Netherlands 



For materials with sufficiently strong uniaxial anisotropy, 
classical elasticity theory admits double-peaked profiles 
with both peak widths and the separation between peaks 
growing linearly as a function of depth . The domain 
of applicability of classical elasticity theory to granu- 
lar materials is not well understood, however, as it of- 
fers no simple way to incorporate noncohesive forces be- 
tween material elements or history dependent frictional 
forces. Several alternative theories for granular stress re- 
sponse have been proposed that make predictions quali- 
tatively different from conventional expectations. Mod- 
els of isostatic materials [l^, and models employing 
"stress-only" consititutive relations [l^l, give rise to hy- 
perbolic differential equations for the stress and predict 
stress propagation along characteristic rays. Similarly, 
the directed force chain network model predicts two dif- 
fusively broadening peaks developing from a single peak 
at shallow depth [2^ . Numerical studies in small isostatic 
or nearly isostatic packings also find evidence of propa- 
gating peaks 0, [l0[- Simulations of weakly disordered 
hexagonal ball-and-spring networks, a common example 
of an elastic material, can display two-peaked stress re- 
sponse when the springs are one-sided [JH^ and uniaxial 
anisotropy is induced by contact breaking. Response in 
the ball-and-spring networks becomes single-peaked as 
friction increases, a result mirrored by a statistical ap- 
proach to hexagonal packings of rigid disks [l^, [l^ . Fi- 
nally, a continuum elasticity theory with a nonanalytic 
stress-strain relation at zero strain has been shown to ac- 
count quantitatively for single-peaked stress response in 
rain-like preparations of granular layers [27j . 

We show here that an elasticity theory incorporating 
both hexagonal anisotropy and near-field microstructure 
effects can account for the experimental observations of 
Geng et al. [H, 0| The theory is phenomenological; it ac- 
counts for the average stresses observed through a com- 
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pilation of many individual response patterns. Our goal 
is to determine whether the ensemble average of effects of 
nonlinearities associated with force chains, contact break- 
ing, and intergrain contact forces can be captured in a 
classical model, and, in particular, to account for the dra- 
matic effects observed in experiments on 2D hexagonally 
close-packed systems. To that end, we develop a non- 
linear continuum elasticity theory applicable to systems 
with hexagonal anisotropy [11] . We find that these effects 
can account for the quantitative discrepancy between the 
Boussincsq solution in 2D (the Flamant solution) for lin- 
ear systems and the experimental data of Refs. [l| and 0| 
for disordered packings of pentagonal grains and hexag- 
onal packings of monodisperse disks. To compare com- 
puted stress fields to the experimental data, we calculate 
the pressure in the material as a function of horizontal 
position at fixed depth. Wc call such a curve a "response 
profile." 



Wc find that induced dipole and quadrupole terms, 
which wc attribute to microstructure effects near the ap- 
plied force, can account for the narrowness of the re- 
sponse profiles in isotropic materials without resorting 
to nonlinear effects. In contrast, the response profiles 
observed in hexagonal packings cannot be fit by the lin- 
ear theory; inclusion of nonlinear terms capable of de- 
scribing hexagonal anisotropy is required. Using a the- 
ory based loosely on a simple triangular lattice of point 
masses connected by springs, but allowing an adjustable 
parameter specifying the degree of hexagonal anisotropy, 
we find reasonable fits to the response profile data. We 
find that for sufficiently strong anisotropy the fitted re- 
sponse profiles correspond to small strains. Thus the 
nonlinear terms are necessary to capture the effects of 
material order, rather than large displacements. This is 
consistent with the experimental observations of Ref . [l[ , 
for which the deformations were small and reversible. 



The paper is organized as follows. In Section [TTl wc 
review well known elements of the theory of nonlinear 
elasticity and the multipole expansion of the stress field. 
In Section IIIIl we develop expressions for the free en- 
ergies of isotropic and several model hexagonal mate- 
rials, including a model in which strong nonlinearities 
arise for small strains. (We use the term "free energy" 
to maintain generality, though in the context of gran- 
ular materials, finite temperature effects are negligible 
and our explicit models make no attempt to include en- 
tropic contributions.) In Section llVl we present a pcr- 
turbative expansion of the response profiles for nonlin- 
ear systems in powers of inverse distance from the point 
of application of the boundary force. In Section |Vl we 
present the response profiles obtained by adjusting the 
monopole, dipole, and quadrupole strengths and the de- 
gree of hexagonal anisotropy. 
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FIG. 1: Stress response in an elastic half-space. Forces must 
vanish everywhere on the free boundary, B, except at the 
origin. The total force transmitted across the surface C is fz. 



II. REVIEW OF ELASTICITY CONCEPTS, 
DEFINITIONS, AND NOTATION 

We first provide a brief review of stress response in 
linear elasticity theory for an isotropic half-plane. We 
then describe the general equations of nonlinear elasticity 
that are solved in Section IIVI for particular forms of the 
free energy. Finally, we review the multipole formalism 
that is later used to model the effects of microstructure 
in the region near the applied force where the continuum 
theory must break down. 

The response of an elastic half-space to a point force 
normal to the boundary, depicted in Fig. [TJ was first 
given by Boussincsq [l7|. A normal force / is applied 
at the origin. In linear elasticity the stress components 
artp and cr^^ vanish on the surface B. The force trans- 
mitted across a surface C enclosing the boundary force 
and with outward normal n must be equal to the force 
applied at the boundary, namely dC z ■ cr ■ fi ~ f . and 

dC X ■ (T ■ fi = 0. We expect that the Boussincsq result 
applies far from the point of forcing, where the stress is 
weak and can be averaged over a large representative vol- 
ume of grains. In this regime, the stress tensor <t is solely 
radially compressive, independent of bulk and shear mod- 
uli, and (in two dimensions) inversely proportional to the 
distance from the point of application 

2/cos0 

cr,.,. = — , ffj-A ~ 0, 0-00 = 0. (1) 

Trr ■ 

Here r and are polar coordinates, being measured 
from the vertical as depicted in Fig. [T] Compressive 
stress is positive. The stress contours are circles passing 
through the origin, where the boundary force is applied. 
This result is a useful approximation to the response in 
a real material far from other boundaries. For linear 
systems, it can be used to calculate the response to an 
arbitrary distribution of force on the boundary. 

Nonlinearities arise from the proper geometric treat- 
ment of finite strains and rotations as well as possible 
anharmonicity in the free energy of the system. In classi- 
cal elasticity, a linear constitutive relation (e.g. Hooke's 
law [2^ ) between stress and strain results from a free en- 
ergy A that is quadratic in the components of the strain 
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tensor. This can be regarded as the first term in a Taylor 
expansion of A about an equiUbrium reference configura- 
tion, and in this paper we include cubic and quartic con- 
tributions to the free energy as well. Unlike the quadratic 
terms, the higher order contributions can distinguish be- 
tween a hexagonally anisotropic system and an isotropic 
one. 

When cubic and higher order powers of the strain 
in A become important, it may also be necessary to 
take into account geometric sources of nonlinearity. Let 
X = {X, Z) be the position of a material element in the 
reference (undeformed) configuration and let x = {x, z) 
be the position of the same material element in the de- 
formed configuration. The displacement field is defined 
as u = X — X and the deformation gradient is defined as 

F = l + Gradu, (2) 

where Grad — [dx ,dz)- To ensure invariancc under over- 
all rotations, one must work with the full Lagrangian 
strain 

r, = \ (F^F - 1) (3) 

rather than just the linearized strain e = (F^"" + F)/2. 
In conventional (linear) elasticity theory, the terms in rj 
nonlinear in u arc neglected and Grad can be replaced 
by grad = {d^,dy,). 

The Cauchy stress cr is the stress measured in experi- 
ments and is a natural function of x. It must satisfy the 
equations of force balance, div cr + pg = 0, and torque 
balance, cr""" = <t, for any deformation. Here div (Div) 
is the divergence with respect to the deformed (unde- 
formed) coordinates. In the context of nonlinear models 
with boundary conditions expressed in terms of forces, 
these equations are more conveniently expressed with re- 
spect to the undeformed coordinates, the nominal stress 
S = JF~^cr, and the reference density po(X) = Jp(x), 
where J = detF. The equations of force and torque 
balance can be rewritten 

DivS + po9 = 0, (4) 
{FSf = FS. (5) 

Defining the thermodynamic tension T via S = TF""", 
the equations are closed by a constitutive relation cou- 
pling T to the Lagrangian strain (and through it the de- 
formation gradient), namely T = Combining these, 
the nominal stress can be written as 

Together, Eqns. ([2][6]) represent a set of equations speci- 
fying the displacements in the system, for a specific ma- 
terial specified by the free energy A, and subject to the 
boundary conditions that stresses vanish on the deformed 
surface (except at the singular point) and the total force 
transmitted through the material is fz. 



Studies of the nonlinear Boussines q pr oblem have fo- 
cused primarily on stability analysis [30l . [3ll . [s^ . Here 
we emphasize the form of the stress response profile and 
restrict our attention to two-dimensional isotropic and 
hexagonally anisotropic systems. As will be described 
below, the stress response can be developed in an expan- 
sion in inverse powers of the distance from the boundary 
force, reminiscent of a multipole expansion of an electro- 
magnetic field. 

The stress response of a hexagonal packing in Ref. [H 
(reproduced in Figs. FTl fTU)) displays a rich structure, de- 
veloping new peaks with increasing depth that gradually 
broaden and fade. It is apparent that Eq. ^ can never 
recreate such a response profile, as there is no length scale 
over which the response develops. However, it is possi- 
ble to create two- (or more) peaked response in isotropic 
linear elasticity. All that is necessary is the application 
of more than one force at the boundary. Two boundary 
forces oriented at ±7r/6 to the normal, for example, will 
produce a two-peaked stress response at shallow depths, 
as shown in Fig. For depths much greater than the 
distance between the two forces, the response approaches 
that of a single normal force equal to the sum of the nor- 
mal components of the two boundary forces. 

At distances larger than the separation between the 
points of appfication of the force, the stress field in Fig. [2^ 
can be closely approximated by a multipole expansion. In 
a granular material, the local arrangement of grains in re- 
gions where strains are large will induce deviations from 
the continuum theory, and in the Boussinesq geometry 
the far field effects of these deviations can be approxi- 
mated by placing a series of multipolar forcing terms at 
the origin. Thus, although the physical force applied by 
Geng et al., for example, was a single, sharply localized, 
normal force, we include in our continuum theory param- 
eters specifying dipole, quadrupole, and perhaps higher 
order multipole forcing strengths to account for the ef- 
fect of microstructure. If the applied force is spread over 
enough grains that the continuum solution predicts only 
small strains everywhere, then the multipole contribu- 
tions can be explicitly computed within the continuum 
theory. If, on the other hand, the force is applied to a 
single grain and represented delta-function in the 
continuum theory, the theory will predict large strains 
near the origin and microstructure effects must be taken 
into account either phenomenologically, as we do here, or 
through a more detailed model of the microstructure in 
the vicinity of the applied force. We conjecture that the 
size of this region near the origin scales with the "iso- 
staticity length scale" discussed in Refs. [l^ and [l^ . 

The first several multipole forces and corresponding 
pressure profiles, are depicted in Fig. [2]3-g. A multipole 
force with stresses that decay as 1 /r" can be constructed 
from n evenly spaced compressive or shearing boundary 
forces having alternating directions and magnitudes in 
proportion to the n"^ row of Pascal's Triangle. The in- 
tegral dx x"'~^ f{x) is the lowest order nonvanishing 
moment of the boundary force distribution f{x). 
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(a) 



ization to arbitrary n: 




FIG. 2: (a) Contour plot of pressure for two point forces 
of equal magnitude located at ±e/2 and oriented at ±7r/6 
from the surface normal. Distances are in units of e. The 
response is two-peaked for shallow depths, transitioning to 
the circular contours of Orr for a single normal force at the 
origin. Monopole (b,e), dipole (c,f), and quadrupole (d,g) 
boundary forcings, along with contours of the corresponding 
pressures. 



The form of the far-field stress response to multipole 
forcing in linear elastieity can be developed by consider- 
ing the Airy stress function x such that cr^r = drxl'"' + 
d^^x/'r'^i 0-1-0 = = -dr{d(i,x/r), and a^^ = drrX- 
The Airy stress function is biharmonic: 



rd0(rsin(/))"-i (p • t*" 
rd(j) (r sill (j))^^^ (^q ■ cr 



(n) 



= 



— kc 



(9) 



where p = x (z) and q = z {x) for odd (even) powers 
of n. k and a carry the units of stress and length, re- 
spectively, the applied force is / = ka. Subject to this 
normalization, the dipole stress tensor cr(^) is 



r(2) - 



.(2) 



.(2) 



2 COS 20 

2 sin 20 



0, 



(10) 



and the quadrupole stress tensor cr^^) jg 



(3) _ 
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5ka^ 

5- cos 30 

7rr3 


3fca3 
5- cos 
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3fca3 

5- sm ( 
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(3) _ 

4,4, - 


— T cos Sep 4 
7rr3 


3ka^ 

^ cos 4. 

7rr3 



A Ax 0. 



(7) 



(11) 



Contours of the associated pressures p*^"' = (l/2)Trcr("' 
and sample boundary forces which produce them are 
shown in Fig. ^jp-d. 

The higher order multipole terms decay more quickly 
than the monopole term, so at asymptotically large depth 
in a material in which both monopole and higher order 
terms are present, the response is indistinguishable from 
the Boussinesq solution. Closer to the point of applica- 
tion, the induced multipole terms contribute more com- 
plex structure to the response. The distance over which 
this structure is observable depends on the material prop- 
erties through the elastic coefficients and increases with 
the strength of the applied force /. 



Assuming x has the form 



X(r,0)=r2^-X(")(0) 



(8) 



and solving for x yields a series of corresponding 
tensors cr^"^ (It is convenient to restrict ourselves to 
transversely symmetric multipole terms, such as those 
in Fig. [IJj-d, so that there is only one corresponding 
stress tensor for each value of n.) cr^^^ corresponds to 

the monopole of Eq. ([1]). For each cr^"\ a^^^J and cr,',^' 
must vanish on the surface except at the origin. For the 
surface C in Fig. [T] we generalize the monopole normal- 



Ill. MODEL FREE ENERGIES 

Here we develop expressions for the elastic free energy 
of several model systems having hexagonal symmetry. 
These will be needed to construct constitutive relations 
relating stress and strain. 

A. Symmetry considerations 

To linear order the elastic energy is quadratic in the 
strain components: 



A = -y^ijki rjij Vki- 



(12) 
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A is a fourth order tensor of rank two, and its compo- 
nents are the elastic coefficients of the material. For 
an isotropic material the free energy must be invariant 
for rotations of rj through arbitrary angle. Therefore A 
can depend only on scalar functions of the strain ten- 
sor components. In two dimensions, the strain tensor 
has two eigenvalues or principal invariants. All other 
scalar invariants, including the independent invariants 
Ii = Trrj = r]ii and I2 = Trrj^ = iVij)"^ (summation 
implied), can be expressed in terms of the principal in- 
variants (ssj or, equivalently, in terms of Ii and l2- The 
free energy of an isotropic linear elastic material can be 
expressed in terms of combinations of Ii and I2 that are 
quadratic in the strain components. 

A=^Xlf+fil2 (13) 




-0.3 -0.2 -0.1 0.1 0.2 0.3 



FIG. 3: (a) A ball-and-spring network with hexagonal symme- 
try and springs oriented horizontally. Even for a linear force 
law, the free energy has terms of cubic and higher order in the 
strains when the equilibrium length of the springs is nonzero, 
(b) Free energy as a function of strain for a unit cell of the 
ball-and-spring network in (a), (solid black) Vertical uniaxial 
compression: rj = rjzz with riyy = 0. (dashed black) 77 = "q^cx 
with Horizontal uniaxial compression: rj^z — 0. (dashed gray) 
Linear elastic approximation for both cases. 77 < corre- 
sponds to compression. 



where A and fi are the Lame coefficients. The reasoning 
generalizes to higher orders. At each order, there will 
be as many elastic coefficients as there are independent 
combinations of Ii and I2. To quartic order in the strains, 
we have 
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UJlIl + LJ2I1I2 



(14) 



We refer to the w's and the Si's as third and fourth order 
elastic coefficients, respectively. 

To construct the free energy of a hexagonal material, 
it is useful to consider a change of coordinates 



material is, to quartic order, 
A = 



^ (^2Ai7?^577Cc + 4A2?7|^^ 



1 



where 7]^^ ^ 77^3, - 77^^ -I- 2177^;^, 77,^^ = ri^x ~ Vzz - ^ir]^^, 
and ?75^ = T]a:x + rjzz- 

For simplicity, we have assumed that terms involving 
gradients of the strains are negligible [13, [s^ . 



B. Hexagonal ball-and-spring network 



^ = X + iz 

C — X — iz, (15) 

as suggested in Ref. [2^. For a rotation of 7r/3 about 
{z X x) these coordinates transform as ^ ^ ^e'^'^^ and 
C Ce-""'/^. The free energy of an elastic material must 
be invariant under such a rotation, which implies that a 
component of the tensor A can be nonzero if and only 
if it too is invariant. For example, the quadratic coef- 
ficient Ajj^^ is nonzero because, under rotation by tt/B, 
A«CC e"/3e'^V3e-'^V3e-'^V3A^^^^ = A^ecc- The only 
other independent nonzero quadratic coefficient is A^^;^^. 
Cubic and higher order coefficients, which arc labeled by 
six or more indices, can also be invariant by having six 
like indices, as in Aj^^^jj. There are three independent 
coefficients at cubic order and four at quartic order. 

The general form of the free energy of a hexagonal 



We now construct the free energy for several specific 
hexagonal materials, taking the point-mass-and-spring 
network of Fig. [3^ starting point. The elastic co- 
efficients are determined by calculating the free energy 
under a homogeneous strain and comparing to Eq. p6|) . 
The springs are taken to have an equilibrium length £ 
and to obey Hooke's law: for a spring with one end at xi 
and the other at X2 the force is 

/ = -fc(V(x2-xi)-(x2-xi)-^), (17) 

where k is the spring constant. We take the springs to 
be at their equilibrium lengths in the undeformed system: 
£ ^ a, the lattice constant. 

Consider the homogeneous strain 
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which stretches the coordinates to a; = y/1 + 2rixx X and 
z = y/1 + 2rjzz Z. The free energy per unit (undeformed) 
volume of a hexagonal ball-and-spring network with one 
third of the springs oriented along the x direction under 
this stretch is 



A 



9 3 



11 3 



9 



135 
128 



Vxx 



215 
^ 128 
45 3 



45 2 2 

1" Q^VxxVzz 



5 

+ 32^^^ 



aVxxVz 



(19) 



The presence of cubic and higher order terms in the 
free energy is due to the nonzero spring equilibrium 
length. The free energy for a constrained axial com- 
pression/extension in the x and z directions is plotted in 
Fig. [21 The corrections to the quadratic expression stiffen 
the system under compression and soften it slightly under 
small extensions. 

Comparing Eqs. pB]) and and equating like coef- 
ficients of r]xx and rj^z wc find 

\/3 7 



Ai — A2 



Ai = -A2 = -A3 = A#fc 



32 



(20) 



Li — L2 



-Li 



15v2l 
128 



A similar calculation for a material in which one third 
of the springs arc oriented vertically, corresponding to 
a reference configuration rotated by 90° from the one 
shown in Fig. [31 yields 



Ai = A2 = 

-Ai = -A2 = -A3 = ^k 
Li = L2 = is = L4 = ^i^ k 



(21) 



1. The a-material 

Goldenberg and Goldhirsch [1, [2^ find two-peaked 
stress response in numerical simulations of a hexago- 
nal lattice of springs when the springs are allowed to 
break under tensile loading. Contact-breaking explicitly 
breaks our assumption of local hexagonal anisotropy in 
any particular sample. In the context of an ensemble av- 
erage, however, the material description retains hexago- 
nal symmetry and the effects of contact breaking are cap- 
tured phenomenologically by considering material made 
of springs with a force law that softens under extension. 



f = -k (V(X2 -Xi) • (X2 



- xi I - a 



+k- 



X2 - xi) • (x2 - xi) - a .(22) 





a = 8 


compressive 


= 4 

tensile « - U 



-0.12 -0.06 0.06 0.12 
6x/a 

FIG. 4: The force law of Eq. dH]) for = 1 and a = . . . 8 



For a > the springs soften under tension and stiffen 
under compression, as shown in Fig. [H In the horizontal 
orientation the clastic constants from Eq. ([^0)) are shifted 
according to 



Ai — A2 



V3l 



Ai = -A2 = -A3 = ^fc-f^fc 



a 3V3 1 



(23) 



Li — L2 — L3 — — L4 — 



_ 15n/3 



h, _ a 9\/3 jL 
128 a 64 



2. The P -material 

In the spirit of phenomenological modeling, all of the 
elastic constants consistent with hexagonal symmetry 
should be considered to be parameters to be determined 
by experiment. To probe the importance of hexagonal 
anisotropy, we consider a model in which all elastic con- 
stants but one arc fixed and define a parameter (3 corre- 
sponding to the strength of the anisotropy. Note that the 
elastic constants for the two orientations of the hexagonal 
ball-and-spring network considered above can be rewrit- 
ten as 



Ai = A2 = ^k 
A2 = A3 = -^fc 



Ai =/3^A: 



(24) 



15\/3 ; 
128 



128 



The case /3 = 1 gives the network with horizontal springs; 
/3 = — 1 gives the network with vertical springs; and /3 — 
gives an isotropic system. Linear response for elastic 
materials with other ansisotropies is treated in Ref. p^ . 



IV. METHOD 

We wish to obtain corrections to the linear elastic re- 
sult for a material with hexagonal symmetry. For later 
convenience we write / = Qmka, where Qm is dimension- 
less, k has units of a spring constant, and a is a lattice 
constant with units of length. We expand the stress in 
successive inverse powers of the radial coordinate, and 
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refer to the terms in the expansion as the dipole correc- 
tion, quadrupole correction, and so forth. For simpHcity 
and clarity, we present here the calculation corresponding 
to the free energy of Eq. (jl6p with coefficients given in 
Eq. 120]) in detail. General equations for arbitrary elastic 
coefficients are exceedingly long and unilluniinating. 

We solve for the the displacements uu{R,^) and 
$), from which the stress tensor can be recon- 
structed. Capitalized coordinates are used as we are now 
careful to distinguish between the deformed and unde- 
formed states. After the deformation, the point X is at 
X = X + uii{R, ^))R + 4))l>. To linear order and 

for the ball-and-spring network described in Eq. (PD|) the 
displacements are 



find 

11 - 13cos2$ - 3cos4<I> - 9cos6$ - 6cos8$ 
= ^v'^ - 27vi - 27w'o + 9w[ 

+ Q<-27u;'i^ln(i?/i?o); 

-5 sin 2$ + sin 4$ + 3 sin 6$ + 2 sin 8$ 
9 

= 3^0 -I- 3v[ + -Wq — 3wi 

+ (3v[ + ^w'A\niR/Ro). (28) 



u^°\r,^) 



"^^'"^ ( cos$ln(i?/i?o) + -$sin$ 

TT V 3 



ViQr 



sin$ln(i?/i?o) 



2 1 

-- sin$ + -^cos'l' ) . (25) 

3 3 



The parameter Rq requires comment. Because the ma- 
terial is semi-infinite in extent, it is free to undergo an 
arbitrary rigid-body translation in the Z-dircction under 
the influence of a normal boundary force. Thus the point 
along the .Z-axis at which the deformation u is zero may 
be chosen arbitrarily. Rq parameterizes this variation. 
Note that the nominal stress, which in the linear theory 
is equivalent to cr in Eq. ([1]), is independent of Rq. 



To find the dipole correction, we take un 



,(0) 



,(1) 



and M$ 



,(0) 



,(1) 



and assume a correction of the form 



u<^^\R,<f) 
4'^(i?,$) 



a — ha — — m(n/Ro) 



R 

jWo(£) 
R 



a 



R 
R 



ln(i?/i?o). (26) 



The deformation gradient F in polar coordinates is 



I + dfiuji {d,s,uji — u^) / R 
Oru^ 1 + (3<i)U$ + uji) /R 



(27) 



Through Eqs.[3]and[6]thc nominal stress can be written 
entirely in terms of the displacements, and through them 
in terms of the four unknown functions, up, wi, uig, and 

Wi. 

Substituting the linear Boussinesq solution of Ea. [25l in 
Eq. mi evaluating Eq. and requiring the coefficient 
of 1/R^ to vanish yields conditions on the w's and w's. 
(Terms of smaller order in l/R vanish identically.) We 



For the moment, we neglect terms of higher order in l/R. 
The source terms on the left-hand side in Eq. ((28l) are 
generated by the linear solution. Requiring coefficients 
of In R to vanish independently gives four second-order 
ordinary differential equations for the four unknown func- 
tions. 

The conditions that normal and shear forces vanish ev- 
erywhere on the deformed boundary except at the point 
of application of the external force can be written 



5M(i?^0,$ = ±7r/2) = 
5<i,cE,(i? 7^ 0,$ = ±7r/2) = 0. 



(29) 



Both the 5*$^ and 5*$$ components of stress have terms 
proportional to Ini?. When we require these terms to 
vanish independently of all other terms, Eq. (p9)) rep- 
resents eight constraints. The nominal stress must also 
satisfy force-transmission conditions 



x-S' -ndC = 



(30) 



where C is any surface enclosing the origin (see 
e.g. Fig. [T|) and n is the unit normal to C. Eq. ([50]) 
is satisfied by the linear elastic solution, and all solutions 
to Eq. ((28)) subject to Eq. ((29)) contribute zero under the 
integration, so this provides no additional constraint on 
the system. 

The eight constraints of Eq. ([29]) fix only seven of the 
eight integration constants. The eighth integration con- 
stant, which we label Qd, multiplies terms identical to 
those contributed in linear elasticity by a horizontally 
oriented dipole forcing such as that depicted in Fig. ^ 
and given in Eq. ((TH)) . Qd is fixed by demanding that a 
variation of the parameter Rq produce only a rigid body 
translation of the material. The integration constants 
determined in this way produce a nominal stress S inde- 
pendent of Rq, as must be the case. 
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FIG. 5: Two imagined scenarios in which a point force induces 
a dipole. Regions of overlap indicate a compressive contact, 
(a) The disks in the second layer slide outward, e.g. for low 
friction, (b) Alternatively the disks might roll inward, to- 
wards the line of force, e.g. due to greater friction between 
grains. This would select a dipole term in the stress response 
with opposite sign from the case depicted in (a). Thus, the 
details of the near field response depend on the mechanics of 
the discrete system. 



The solution of Eq. [28] consistent with Eq. ([291) is 
"5 7 



uo($) 



6 ' 3 
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■ cos 2$ 
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■ COS 4$ + - cos 6$ 
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— $ + - sin2$ 
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/ 47rQrf 
V 73 

2tt' 
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+ ^ sin 4$ -[ 






' 12 
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11 


+ — sin8$ - 




36 


18* 


sin 2$ 





27r 



2 , sin2$ln(i?/i?o). 



(31) 



For the choice Rq = a, we find the induced dipole 
coefficient Qd = 0, and for the sequel we fix i?o to have 
this value. The same choice of i?o also yields the induced 
quadrupole coefficient Qq = Q below. As discussed above, 
rather than set them to zero, we leave these terms in the 
displacements, and correspondingly the stresses, as free 
parameters to account for the influence of microstructure 
on the response. They are weighted so that Qd ~ 1 and 
Qq = 1 correspond to the stresses of Eqs. ([T0[) and pT|) . 



We repeat the process described above to develop 
quadrupole corrections to the stress response. The dis- 
placements arc assumed to have the form ur{R, $) = 



(7?,$) 



,(2) 



and w<i,(i?,$) = 
(i?, $) where the second or- 



der corrections have the form 



,(2) 



,(2; 



v^m Vim 

i?2 



i?2 

Worn 



i?2 

In {R/Rq) 



In (i?/i?o) 



i?2 

i?2 



i?2 

In {R/Rq 



In {R/Ro) 



(32) 



The details of the calculation are omitted, as they are 
conceptually similar to the dipole calculation but involve 
much longer expressions. Defining c„ = cosn<f>, s„ = 
sinn<i>, and L = In {R/ Rq), the pressure is 
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(33) 



where B2 = QlM^ j Yl^/Z j t^"^ , B3 = Qf„ka^/36Tr^, and 
B's = 4Q^fca3/3%/37r2. 

We will find below that the /3-niaterial best describes 
the data of Ref . P| . In this case the pressure of Eq. ([33]) 
gains a number of additional terms involving f3. These 
terms are given in the Appendix. 



V. RESULTS 

Given the expressions derived above for the pressure, 
we perform numerical fits to the data from Geng et al. [l[ . 
There are four fitting parameters for the ball-and-spring 
material: the monopole coefficient Qm, the dipole coef- 
ficient Qd, the quadrupole coefficient Qq, and the spring 
constant k. We take the lattice constant to be the disk 
diameter: a = 0.8 cm. The three multipole coefficients 
have been defined to be dimensionless. We set = a 
so that Qd and Qq would be zero in a theory with no 
microstruturc correction. In two dimensions the units of 
stress are the same as the units of the spring constant k. 
Thus k sets the overall scale for the stress. For theoretical 
purposes, k could be scaled to unity; in our fits it serves 
merely to match the units of stress in the experimental 
data. 

We attempt to fit experimental measurements on 
pentagonal-grain packings by varying Qm, Qd and Qq 
in the isotropic theory. To explain the experimental data 
on hexagonal disk packings, wc attempt fits based on 
the ball-and-spring network, the a-material, and the /3- 
material. 

Wc regard the response profiles presented in the fol- 
lowing section, particularly Figs. [S] and [51 as a proof of 
principle: average response in experiments of the sort 
performed in Ref. [l| is consistent with an elastic con- 
tinuum approach when microstructure and material or- 
der are properly incorporated. The results we present 
arc phcnomonological in that we have obtained elastic 
coefficients and multipole strengths by fitting to data. 



We expect that the elastic coefficients we fit are material 
properties in the sense that they could be determined by 
experiment or simulation in another geometry (e.g. a uni- 
form shear or compression) , then used in our calculations 
for point response. 



A. Fitting to pressure 

The photoelastic measurements of Geng et al. asso- 
ciate a scalar quantity with each point in space. The 
measurement technique extracts no directional informa- 
tion, so the relevant theoretical prediction to compare to 
experiment is the local pressure p = (l/2)Tr cr [36[. 

The data of Ref. [ij are averaged over many experi- 
mental realizations; the average hydrostatic head is also 
subtracted. The hydrostatic contribution to the stress 
is largest at depth where, as seen below, the linear 
(monopole) response dominates. Therefore, although the 
elasticity theory is nonlinear and superposition docs not 
strictly hold, we expect the incurred error from differ- 
encing to be small. We note also that our fits necessarily 
produce regions of small tensile stress near the surface. 
Removal of all tensile stresses from the solution would 
require treating the nonlinearity associated with contact 
breaking to all orders in the nonlinear elasticity theory. 
In the present context, such regions should be taken only 
as indicating that contacts are likely to break. 

Fitting to the Cauchy pressure p, which is a nat- 
ural function of the deformed coordinates x, presents 
a difficulty. Namely, our caluclations yield a relation 
X = X + u(X) that is not invcrtiblc. Although in princi- 
ple cr is known for all points in the deformed material, we 
can still only reference those points by their undeformed 
positions. That is, we have calculated p(x(X)). Thus for 
the purposes of fitting, we neglect the difference betwen 
X and X. In the experiment, the forcing was restricted to 
strengths for which the strains were small; there were no 
large-scale rearrangements. This suggests that replacing 
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the deformed coordinates with the undeformed coordi- 
nates will introduce only small errors. Of course, if the 
strains are small, it is reasonable to ask whether nonlin- 
ear elasticity is really needed or helpful. A discussion of 
this point is provided in Section IVll below. 

To facilitate comparison between various materials, we 
restrict our consideration to boundary forces / = Qmka 
with Qm = 1- We have found that similar response pro- 
files can be obtained for 0.25 < Qm < 2, and all best-fit 
values for Qm lie in this range. The force f = ka is that 
required to compress one Hookean spring through one 
lattice constant. 

Rather than compare pressure directly to the data of 
Rcf. [l|, we scale each data point by its depth Z and 
fit to ZP{X,Z) for two depths: Z = 2.7 cm and 3.60 
cm (recall that the grain diameter is 0.80 cm). Scaling 
by Z compensates for the decay of the response with 
depth. For a reasonable fit, fitting to data at one or 
two shallow depths gives good agreement with all data 
at greater depth. Generally the fitting algorithm returns 
parameters such that agreement with experimental pro- 
files at depths shallower than the shallowest fitting depth 
is poor. For the best model material, however, it is possi- 
ble to achieve reasonable agreement with data at a depth 
of 2.25 cm. 
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depth = 2.70 cm 



B. Pentagonal particles 
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FIG. 6: (color online) (black curves) A fit of Cauchy pressure 
for a spring-like isotropic (/3 = 0) material with free energy 
expanded to quartic order in the strains. The fit parameters 



are Q„ 



1, Qd = 0.5, 



-4.6, and k = 702 and were 



determined by fitting to response in a packing of pentagonal 
particles (gray points) of width 0.8 cm at depths Z — 2.7 cm 
and 3.60 cm. (dashed green curves) Linear elastic multipole 
response with Qm = 1, Qd ~ 0.6, Qq = —4.0, and k = 700, fit 
to the same data, (dotted red curves) Linear elastic monopole 
response with Qm = 1 and k = 1032. 



The nominal pressure of the spring-like isotropic (/3 = 
0) material for Qm ~ Qd ~ 0.5, Qq = —4.6, and 
k = 702 is shown in Fig. [6l Parameters were deter- 
mined by fitting to mean pentagonal particle response 
data. The result is a clear improvement over the fit to 
linear elastic pressure; the nonlinear calculation is able 
to capture the narrowing of the response as Z — > 0. At 
Z = 2.25 cm, shallower than the fitting data, the curve 
has an appropriate width but overshoots the peak. Note 
that there is little reason a priori to assume the elastic 
coefficients we have chosen are the appropriate ones to 
describe this material. 

A multipole expansion 



,ka 



ttR 



■ cos $ ' 



4:Qdka'^ 



cos2$- 



2Qqka^ 
ttR^ 



cos 3$ (34) 



with Qm ^ 1, Qd ^ 0.6, Qq = -4.0, and k = 700 is 
nearly indistinguishable from the full nonlinear expres- 
sion with microstructure correction. This suggests that 
in the disordered packings the deviation from monopole- 
like linear elastic response is a consequence of microstruc- 
ture, not effects captured by the nonlinear theory. 



C. Hexagonal packings 



1. Ball-and-spring fit 



The nominal pressure of the ball-and-spring network 



for Qm = 1, 



9.1, Qq = 36 and k = 112 is shown in 



Fig. [71 Parameters were determined by fitting to mean 
hexagonal packing response data. The pressure has two 
peaks at shallow depths; by Z = 5 cm it has crossed 
over to a single central peak. As expected, the elastic 
prediction improves with depth, as the monopole term, 
which is independent of all elastic coefficients, comes to 
dominate. For depths ^ ^ 3 cm there are clear qualitative 
differences between the fit and the data. The two large 
peaks in the data are wider apart than the prediction and 
they fall off more sharply with horizontal distance from 
the center; moreover, the theoretical prediction fails to 
capture the small central peak in the data. 
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FIG. 7: (black curves)A fit of Cauchy pressure for a ball- 
and-spring network including cubic and quartic corrections 
to the free energy. The fit parameters are Q„i = 1, Qd = 9.1, 
Qq = 32, and fc = 112 and were determined by fitting to 
response in a monodisperse hexagonal packing of disks (gray 
points) of diameter 0.8 cm at depths Z — 2.7 cm and 3.60 cm. 



FIG. 8: (black curves) A fit of Cauchy pressure for the a- 
material including cubic and quartic corrections to the free 
energy. The fit parameters are Qm = 1, Qd ~ 0.9, Qq — —15, 
k — 354, and a = 8.9 and were determined by fitting to 
response in a monodisperse hexagonal packing of disks (gray 
points) of diameter 0.8 cm at depths Z — 2.7 cm and 3.60 cm. 



2. a-material fit 

The nominal pressure of the a-material for Qm = 1, 
Qd = 0.9, Qq = —15.4, k = 354 and a = 8.9 is shown in 
FiglHl The pressure response in the a-material is a better 
fit than that for the ball-and-spring network, as it more 
closely recreates the two-peaked structure from Z w 4 cm 
to 6 cm. It also drops off more sharply in the wings than 
the ball-and-spring response. The central peak, however, 
is still absent. Moreover, a value of a ~ 9 is fairly large 
(see Fig. E]). 



3. 13 -material fit 

The nominal pressure of the /3-material for Qm = 1, 
Qd = 0.6, Qq = -2.0, k = 353 and f3 = 12.4 is shown in 
Fig. [9l Parameters were determined by fitting to mean 
hexagonal response data. The /3-material response does a 



better job of capturing the peaks than both the ball-and- 
spring material and a-material response profiles. Like 
the a-material, the shape of the response peaks of the /3- 
material is narrower and more appropriately positioned 
than that of the ball-and-spring material. The /3-material 
profiles do a better job capturing the small central peak, 
though the required (3 value of 12.4 represents a hexago- 
nal anisotropy that is very strong compared to that of a 
simple ball-and-spring network. 

Fig. [TO] shows the /3-material response without mi- 
crostructure correction {Qm = 1, /3 = 10.8, k = 509) and 
the linear elastic response with induced multipole terms 
of Eq. dMl) (Qm - 1, = 11.4, Qq = 42, fc = 116). Nei- 
ther agrees well with the data. It is necessary to include 
nonlinear as well as microstructure corrections to the lin- 
ear elastic result to obtain good agreement with the mean 
hexagonal response data. This contrasts with the mean 
disordered response data, which can be described with 
a microstructure correction alone. We infer that nonlin- 
ear corrections are needed in the hexagonal system to 
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FIG. 9: (black curves) A fit of Cauchy pressure for the /3- 
material including cubic and quartic corrections to the free 
energy. The fit parameters are Qm = 1, Qd = 0.6, Qq = —2.0, 
k = 353, and /3 = 12.4 and were determined by fitting to 
response in a monodisperse hexagonal packing of disks (gray 
points) of diameter 0.8 cm at depths Z — 2.7 cm and 3.60 cm. 



capture the material anisotropy. 



FIG. 10: (color online) (black curves) A fit of Cauchy pres- 
sure for the /3-material including cubic and quartic corrections 
to the free energy but without multipole corrections for mi- 
crostructure {Qd = = Qm). The fit parameters are Qm = 1, 
k = 509, and /3 = 10.8 and were determined by fitting to re- 
sponse in a monodisperse hexagonal packing of disks (gray 
points) of diameter 0.8 cm at depths Z — 2.7 cm and 3.60 cm. 
(dashed green curves) Linear elastic multipole response with 
Qm = 1, Qd = 11.4, Qq — 43, and k = 116, fit to the same 
data. 



D. Crossover to linear elasticity 

For shallow depths the hexagonal anisotropy of the or- 
dered disk packing is strongly reflected in the functional 
form of its stress response. The dipole and quadrupole 
corrections which shape the response in the near field fall 
off as and 1/i?^, respectively, while the monopole 

response decays as 1/R. Sufficiently deep within the ma- 
terial, the monopole term, which is identical to the linear 
elastic solution, will dominate. Fig. I 111 shows contours of 
the nominal pressure for the /3-material of Fig. [5] in the 
near and far fields. In the first 6 cm of depth the three 
peaks seen in the data are clearly visible. The contours 
of the pressure in linear elasticity are circles, and by a 
depth of 40 cm this form is largely recovered. 



E. Physical pressure and strain 

Having determined fit parameters, it is possible 
to visualize the physical or Cauchy pressure p = 
(l/2)Tr(T(x(X)) and strains in the material. In the un- 
deformed material, each disk sits on a lattice site which 
we label by an index i. Under the deformation the disk 
at Xi moves to x; = -I- u^. We draw a disk of ra- 
dius D = 0.8 cm at x; and shade it in proportion to 
|Xj|pi(xi(Xi)). The first three layers of the packing, for 
which the displacements and pressure are clearly diverg- 
ing near the applied force, are drawn but not shaded. 
Though we do not make any attempt to portray the de- 
formations of the disks themselves, the overlap or sepa- 
ration between disks gives a good sense of the strain in 
the material, and the colormap indicates the local varia- 
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FIG. 11: (color on-line) Pressure contours for the /3-material 
with fit parameters identical to those in Fig. |9] At shallow 
depths the structure is three-peaked, the outer two seeming 
to propagate with depth. At greater depth the crossover to 
monopole response is evident. Regions of tensile stress near 
the surface are plotted in green. 



(b) 




FIG. 12: (a) The deformed /3-material. The first three lay- 
ers of the material are omitted. Disk i with lattice position 
Xi in the undeformed material is shown here centered at x^. 
Each disk is shaded according to Ripi, the physical pressure 
scaled by the (undeformed) distance from the point force; val- 
ues increase from blue through purple and orange to green. 
Pressures are calculated for the case for Qm ~ 1, Qd = 0.6, 
Qq = —2.0, k — 353 and /3 = 12.4. Two-peaked structure 
is apparent, as well as arching in the upper layers. The 
strains are large, (b) The deformed /3-material for Qm = 1/4, 
Qd = 0.12, Qg = -1.2, k = 1415 and 13 = 45.7. 



tion of pressure on the grain scale. The /3-niaterial fit for 
Qm = 1 is shown in Fig. [121 The two-peaked response 
structure is immediately apparent; the smaller third peak 
is more difficult to see, but present for the first few rows. 
There is dilatation near the surface. The disks directly 
below the applied force participate in the formation of 
arches, which is consistent with the appearance of two 
large peaks along the lines (j) = ±7r/6. 



VI. STRAIN MAGNITUDE 

We have demonstrated that hexagonally anisotropic 
nonlinear elastic response can display stress profiles sim- 
ilar to those seen in ordered granular packings, which sug- 
gests that significant deviations from the classical Boussi- 
nesq response can extend to depths of tens of layers. 
However, from Fig. [T2k it is also clear that the atten- 
dant strains are large, creating regions of strains in the 
first two dozen layers that are much larger than those 
observed in the systems studied by Geng et al. This is 
not entirely surprising for the choice Q„i = 1. We note, 
however, that by fixing Qm — 1/4, as in Fig. [T^ . we 
obtain a fit in which strains outside the first three layers 
are reasonably small. Differences from the response pro- 
files in Fig. [9] arc imperceptibly small; plotted on top of 
Fig. [51 the Q,n = 1/4 and Qm = 1 curves would overlap. 
The microstructure corrections are still of order unity, the 
spring constant is four times larger (so that the imposed 
force / is unchanged), and the hexagonal anisotropy is 
increased significantly: f3 = 45.7. Thus in a our simplis- 
tic ball-and-spring-inspired material, the observed pro- 
files can be attributed either to strong nonlinearity due to 
large strain magnitude or to strong hexagonal anisotropy. 

The material constants of Eq. were chosen as a 
minimal hexagonally anisotropic model, rather than de- 
rived from a microscopic model. We speculate that the 
enhancement of the nonlinearity and/or the hexagonal 
anisotropy over the values obtained naturally from sim- 
ple ball-and-spring models may be due to the importance 
of a short length scale S <^ D in the grain-grain inter- 
actions. Such a length scale may be the consequence of, 
e.g., nonhnear grain interactions ("soft shell" grains [13] 
or Hertzian force laws), or inhomogeneous elastic coeffi- 
cients due to microscopic grain irregularities [2^, [sl] , in 
which case small strains may correspond to large defor- 
mations of contacts on the relevant scale 6. Full consid- 
eration of such effects is beyond the scope of the present 
work. 

Considering all the results presented above, we arrive 
at the following picture. The important distinctions be- 
tween 2D disordered and hexagonal granular packings 
are the effects near the applied point force and the ma- 
terial symmetry. Although nonlinearity complicates cal- 
culations considerably, it enters only as a matter of ne- 
cessity in incorporating material order: elasticity can- 
not distinguish isotropic and hexagonally anisotropic ma- 
terials otherwise. The facts that 1) nonlincaritics in 
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the isotropic material provide no notable improvement 
over microstructure corrections alone (see Fig. (5]), and 2) 
hexagonal materials admit reasonable response profiles 
for small strain and strong anisotropy (see Fig. 112b ). mi- 
derscore this point. A large f3 value may be difficult to 
interpret in terms of a microscopic model, but this is 
not surprising given that it represents a combination of 
strong local nonlincarites and an ensemble average over 
microstructures that are known to lead to vastly different 
stress or force chain patterns. 

VII. CONCLUSION 

Our results indicate that continuum elasticity theory 
can provide semi-quantitative explanations of nontrivial 
experimental results on granular materials. For isotropic 
(disordered) materials subject to a point force, it ap- 
pears that nonlincarities are less important than mul- 
tipoles induced at the surface where continuum theory 
breaks down. For hexagonal disk packings, however, the 
anisotropy associated with nonlinear terms in the elas- 
ticity theory is required. We have studied the nonlinear 
theory of response in a hexagonal lattice of point masses 
connected by springs and a phcnomenological free energy 
with an adjustable parameter determining the strength of 
the hexagonal anisotropy. A similar treatment would be 
possible for systems with, e.g. , square or uniaxial sym- 
metry, but the free energy would acquire additional terms 
at all orders. For a particular choice of elastic coefhcients, 
the multiple peaks in the pressure profile at intermedi- 



ate depths and the recovery of the familiar single peak of 
conventional (linear) elasticity theory at large depths are 
well described by the theory. To the extent that theoret- 
ical approaches based on properties of isostatic systems 
predict hyperbolic response profiles (23j , our analysis in- 
dicates that the materials studied in Refs. Q and 
have average coordination numbers that place them in 
the elastic rather than isostatic regime. 
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APPENDIX A: PRESSURE IN THE /3-MATERIAL 

The expression of Eq. ([55]) gives the pressure in a hor- 
izontally oriented ball-and-spring network with elastic 
coefficients given by Eq. ([TB)) . The pressure in the 
material of Eq. (|24p is given by adding an additional term 
pp{r{R, ^),(j){R, $)) to the expression in Eq. ([33]). where 
Pfj is given by 



p;3(r(i?,$),0(i?,$)) 



6C4 — 9C6 

28333 



- 4C8 

205 



139 



35 



-C3 



-C5 



420 2 
2$f -62s3 + 33s5 



Cv H Cn 

3 2-^ 
55 

" 44S7 + —Sg 



L{ -50c3 + 36c7 + 30c9 



V37ri?3 



10418 
9ci ^77^ C3 



735 
5V37r 



15 63 119 10 

ycg + —Cn + ^ci3 + —ci5 + 36$S3 + 24Lc3 



14 



-C3 



28333_7r 
6720 



-C3 



2057r 



32 



-C5 



6v37rc7 H — —C7 



5\/37rc9 + 



, 31 33 11 55 

$ S3 H St; H S7 Sq 

1 4 ^ 8 ^ 2 ^ 24 ^ 



^,25 9 15 

L C3 H C7 -i Co 

I 8 ^ 4 ^ 8 ^ 



16 



Cl 



5209 
5880 



C3 



15 
32 



63 
64 



Cn 



119 
160 



Cl3 



24 



9 3 

Cl5 + -^>S3 -I- -LC3 



(Al) 



[1] J. Geng, D. Howell, E. Longhi, R. P. Behringer, Phys. Rev. Lett. 87, 035506 (2001). 

G. Reydellet, L. Vanel, E. Clement, and S. Luding, 



15 



[2] J. Geng, G. Reydellet, E. Clement, and R. P. Behringer, 

Physica D 182^ 274 (2003). 
[3] D. Serero, G. Reydellet, P. Claudin, E. Clement, and 

D. Levine, Eur. Phys. J. E 6, 169 (2001). 
[4] G. Reydellet and E. Clement, Phys. Rev. Lett. 86, 3308 

(2001). 

[5] N. W. Mueggenburg, H. M. Jaeger, and S. R. Nagel, 

Phys. Rev. E 66, 031304 (2002). 
[6] M. J. Spannuth, N. W. Mueggenburg, H. M. Jaeger, and 

5. R. Nagel, Granular Matter 6, 215 (2004). 

[7] D. Head, A. Tkachenko, and T. Witten, Eur. Phys. J. E 

6, 99 (2001). 

[8] C. Goldenberg and I. Goldhirsch, Granular Matter 6, 87 

(2004) . 

[9] C. Goldenberg and I. Goldhirsch, Nature 435, 188 

(2005) . 

[10] A. Kasahara and H. Nakanishi, Phys. Rev. E 70, 051309 
(2004). 

[11] C. F. Moukarzel, H. Pacheco-Martinez, J. C. Ruiz- 
Suarez, and A. M. Vidales, Granular Matter 6, 61 (2004). 

[12] N. Gland, P. Wang, and H. A. Makse, Eur. Phys. J. E 
20, 179 (2006). 

[13] W. G. EUenbroek, E. Somfai, W. van Saarloos, and 
M. van Hecke, in Powders and Grains 2005, edited 
by R. Garcia-Rojo, H. Herrmann, and S. McNamara 
(A. A. Balkema Publishers, Rotterdam, 2005), pp. 377- 
380. 

[14] W. G. EUenbroek, E. Somfai, M. van Hecke, and W. van 
Saarloos, Phys. Rev. Lett. 97, 258001 (2006). 

[15] S. Ostojic and D. Panja, Europhys. Lett. 71, 70 (2005). 

[16] S. Ostojic and D. Panja, cond-mat/0606349 (2006). 

[17] J. Boussinesq, Application des Potentiels d I'Etude 
de I'Equilibre et du Mouvement des Solides Elastiques 
(Gauthier-Villars, Paris, 1885). 

[18] M. Otto, J.-P. Bouchaud, P. Claudin, and J. E. S. Soco- 
lar, Phys. Rev. E 67, 031302 (2003). 

[19] M. Wyart, S. R. Nagel, and T. A. Witten, Euro- 
phys. Lett. 72, 486 (2005). 



[20] R. C. Ball and R. Blumenfeld, Phys. Rev. Lett. 88, 
115505 (2002). 

[21] A. Tordesillas, S. D. C. Walsh, and B. S. Gardiner, BIT 

Numerical Mathematics 44, 539 (2004). 
[22] A. V. Tkachenko and T. A. Witten, Phys. Rev. E 60, 

687 (1999). 

[23] R. Blumenfeld, Phys. Rev. Lett. 93, 108301 (2004). 

[24] J.-P. Bouchaud, P. Claudin, M. Gates, and J. Wittmer, in 
Physics of Dry Granular Media, edited by H. J. Hermann, 
J. P. Hovi, and S. Luding (Kluwer Academic, Boston, 
1997). 

[25] J. E. S. Socolar, D. Schaeffer, and P. Claudin, 

Eur. Phys. J. E 7, 353 (2002). 
[26] C. Goldenberg and I. Goldhirsch, Phys. Rev. Lett. 89, 

084302 (2002). 

[27] K. Brauer, M. Pfitzner, D. O. Krimer, M. Mayer, 

Y. Jiang, and M. Liu, Phys. Rev. E 74, 061311 (2006). 
[28] R. W. Ogden, Non-Linear Elastic Deformations (Dover 

Publications, Mineola and New York, 1984). 
[29] L. D. Landau and E. M. Lifshitz, Theory of Elasticity 

(Butterworth-Heineman, Oxford, 1997). 
[30] J. G. Simmonds and P. G. Warne, J. of Elasticity 34, 69 

(1994). 

[31] E. T. Coon, D. P. Warne, and P. G. Warne, J. of Elas- 
ticity 75, 197 (2004). 

[32] M. R. Lee, D. P. Warne, and P. G. Warne, Mathematics 
and Mechanics of Solids 9, 97 (2004). 

[33] A. J. M. Spencer, Gontinuumm Mechanics (Dover Pub- 
lications, Mineola and New York, 1980). 

[34] A. S. J. Suiker, R. de Borst, and C. S. Chang, Acta Me- 
chanica 149, 161 (2001). 

[35] S. D. C. Walsh and A. Tordesillas, Granular Matter 6, 
27 (2004). 

[36] J. Geng, Ph.D. thesis, Duke University, Durham, NC, 
USA (2003). 

[37] P. G. de Gennes, Rev. Mod. Phys. 71, S374 (1999). 
[38] B. A. DiDonna and T. C. Lubensky, Phys. Rev.E 72, 
066619 (2005). 



